Astronomy & Astrophysics manuscript no. article 


© ESO 2008 


5th February 2008 





Acoustic oscillations of rapidly rotating polytropic stars 
II. Effects of the Coriolis and centrifugal accelerations 

D. Reese, F. Lignieres, and M. Rieutord 

Laboratoire d'Astrophysique de Toulouse et Tarbes - UMR 5572 - Universite Paul Sabatier Toulouse 3 - Observatoire Midi- 
Pyrenees, 14 avenue E. Belin, 31400 Toulouse, France 

Received March 24, 2006; accepted Mai 12, 2006 

Abstract 



o 
o 

<N 
<D ■ 

\o • 

(N 
> 

^ Context. With the launch of space missions devoted to asteroseismology (like COROT), the scientific community will soon have 
ly-^ , accurate measurements of pulsation frequencies in many rapidly rotating stars. 

■ Aims. The present work focuses on the effects of rotation on pulsations of rapidly rotating stars when both the Coriolis and 
f"^) centrifugal accelerations require a non-perturbative treatment. 

\& . Methods. We develop a 2-dimensional spectral numerical approach which allows us to compute acoustic modes in centrifugally 
distorted polytropes including the full influence of the Coriolis force. This method is validated through comparisons with 
previous studies, and the results are shown to be highly accurate. 

Results. In the frequency range considered and with COROT's accuracy, we establish a domain of validity for perturbative 
methods, thus showing the need for complete calculations beyond v ■ sin i = 50km.s _1 for a R = 2.3 R Q , M = 1.9 M polytropic 
star. Furthermore, it is shown that the main differences between complete and perturbative calculations come essentially from 
the centrifugal distortion. 
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03 ' 1. Introduction 

The study of rapidly rotating stars is a field in which there are many unresolved questions. The structure of these 
stars, the rotation profile, the angular momentum transport and many other aspects are not well understood. In order 
to answer some of these questions, many different theoretical and observational methods have been developed over the 
years. For instance, interferometry is starting to give clues as to the shape of these stars and effects such as gravitational 
darkening (e.g. Domiciano de Souza et al., 2003, 2005, Peterson et al., 2006). On the theoretical side, there exists a 
number of numerical models which are progressively becoming more realistic (e.g. Roxburgh, 2004, Jackson et al., 
2005, Rieutord et al., 2005, Rieutord, 2006). These models can then be supplemented with asteroseismology which 
relates the internal structure to observable stellar pulsations. In order to fully exploit these pulsations, it is necessary 
to accurately quantify how they are affected by rotation. In the present work, we will show how this can be done for 
acoustic pulsations in uniformly rotating polytropic stellar models. 

Rotation has several effects on stars and their pulsations. These result from the apparition of two inertial forces, 
namely the centrifugal and the Coriolis forces. The centrifugal force distorts the shape of the star and modifies its 
equilibrium structure. The Coriolis force intervenes directly in the oscillatory motions. Neither of these effects respect 
spherical geometry, which means that the radial coordinate r and the colatitude 9 are no longer separable. As a result, 
pulsation modes cannot be described by a single spherical harmonic as was the case for non-rotating stars. In order to 
tackle this problem, two basic approaches have been developed. The first one is the perturbative approach and applies 
to small rotation rates. In this approach, both the equilibrium structure and the pulsation mode are the sum of a 
spherical solution (or a single spherical harmonic), a perturbation, and a remainder which is neglected. The second 
approach consists in solving directly the 2-dimensional eigenvalue problem fully including the effects of rotation. 
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Historically, the perturbative method has been applied to first, second and third order. Previous studies include 
Saio (1981), Gough & Thompson (1990), and Dzicmbowski & Goode (1992) for second order methods and Soufi et al. 
(1998) and Karami et al. (2005) for third order methods. These have been applied to polytropic models (Saio, 1981) 
and then to more realistic models. There have also been some studies based on the non-perturbative approach. Most 
non-perturbative calculations have focused on the stability of neutron stars, r- modes and f-modes rather than on 
p-modes. Some exceptions are Clement (1981, 1984, 1986, 1989, 1998), Yoshida & Eriguchi (2001) and Espinosa et al. 
(2004). 

The present work aims at accurately taking into account the effects of rotation on stellar acoustic pulsations, 
so as to be able to deduce asteroseismological information from rapidly rotating stars. Previous results are either 
inaccurate or not valid for high enough rotation rates. In order to achieve a sufficient degree of precision, we used 
numerical methods which have already proved to be highly accurate for other similar problems. The present method 
is a further development of the numerical method of Lignicrcs et al. (2006b) and Lignieres et al. (2006a) (hereafter 
Paper I) who used a spectral method (Canuto et al., 1988) with a surface-fitting spheroidal coordinate system based 
on Bonazzola et al. (1998). The spectral method itself has already been used for calculating inertial waves in spherical 
shells (Rieutord & Valdettaro, 1997) and gravito-inertial modes in a 1.5 M Q ZAMS star (Dintrans & Rieutord, 2000) 
both of which involve the non-perturbative effects of the Coriolis force. Achieving high precision is of great importance 
for interpreting present and future measurements of stellar pulsations. Furthermore, it provides a means to establish 
the domain of validity of perturbative methods. Finally, this work can then be used as a reference to validate future 
methods. 

The organisation of the paper is as follows: in the next section, the numerical method is described in detail. This 
section is followed by a series of comparisons and tests which establish the accuracy of the results. We then proceed 
to discuss perturbative methods and their validity. A conclusion and outlooks follow. 

2. Formalism 

The calculation of oscillation modes of rotating polytropes takes place in two steps. Firstly, an equilibrium model must 
be determined. Secondly, this model needs to be perturbed so as to give the eigenoscillations. 

2.1. Equilibrium model 

The equilibrium model is a self-gravitating uniformly rotating polytrope described by the following equations in the 
rotating frame: 

Po = Kpl, (1) 

= -VP - Po V (* c - ±fi 2 a 2 ) , (2) 

A* = 4nG Po , (3) 

where P is the pressure, p the density, K the polytropic constant, 7 the polytropic exponent, ty the gravitational 
potential, s the distance to the rotation axis and G the gravitational constant. One can also introduce the polytropic 
index N = 1/(7 — 1) and a (pseudo-)enthalpy h = J dP/p = (1 + N)P /p . The pressure and density profiles are 
then proportional to powers of this enthalpy: P oc h N+1 , and p Q oc h N . A number of non-dimensional parameters also 
intervene and characterise the polytropic model: 

AnGpcR 2 nR cq Pc i? pol 

A= ft* = -7r^' a =7T' £ = 1- , (4) 

h c y/h c {p) -fteq 

where quantities with the subscript "c" denote the equilibrium value at the centre of the polytrope, i? eq and R po \ are 
the equatorial and polar radii, resp., and (p) = 3M /AttR^ a pseudo-mean density. The method used to compute the 
equilibrium model is described in Paper I. 

2.2. Perturbation equations 

We calculate adiabatic, inviscid oscillation modes using Eulerian perturbations to the equilibrium quantities 1 . The 
linearised equations in the rotating frame read: 

d t p = V • ( Po v) , (5) 



1 The term "perturbation", which means a small departure from equilibrium in this context, is not to be confused with 
perturbation from the perturbative method, where it means a small departure from the spherical case. 
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Po d t v = - Vp + pg - /9 V* - 2 Po n x v, (6) 

n N 2 r 2 

dtp-cZdtp = ^rfvg , (7) 

= A* - AirGp, (8) 

where quantities with the subscript "o" denote equilibrium quantities and those without any subscript Eulerian 

perturbations. g is the effective gravity, c is the speed of sound, T\ the adiabatic exponent and N Q the Brunt- 
Vaisala, frequency. These are given by the following formulas: 

g = - v U a - invV (9) 



2 

c 2 = T.Po/po, (10) 

r - ■ (SSL 

It is worth noting that we have used the fluid's barotropicity in the definition of N a . 

We can then put these equations in non-dimensional form using the following transformations: 



t = T r t, p = p c p, r = R cq r, g = g r g, v = V r v, 

p = P r p, fl = UJrQ Co = V r C Q , N = UJ r N_o, 

where: 



(13) 



W r =T- 1 = (*KG Pc ) 1 ' 2 , V r = ^, g r =^= A -KGR e(l p c , Pr = ^. (14) 

It is important to note that tt± and £2 correspond to two different dimcnsionless expressions of the rotation rate. In 
order to go from one expression to the other, one can use the following formula: 

n*=£h/A, (15) 
where A is given by Eq. (4). 

If wc assume a time dependence of the form exp(Ai), the following generalised eigenvalue problem is obtained (we 
have dropped the underlined notation): 

Xp = -v ■ Vp - p V ■ v, (16) 

Xp v = - Vp + pg Q - p V* - 2p a Qe z x v, (17) 
n N 2 r 2 

\p-\c 2 o p = P^LJL V .g , (18) 

\\9o\r 

= A* - p. (19) 
2.3. Change of variables 

In order to have solutions with a good numerical behaviour on the surface of the star, we use the following variables: 

n= i^' 6= F^' (20) 

where H = h /h c is a non-dimensional form of the enthalpy. These choices result from an analysis of the behaviour 
of the solution near the surface, based on a "generalised" Frobcnius study of the system of equations. Although not 
fully proved, this study gives the correct results in the spherical case (see Appendix A). It also leads to the following 
boundary condition on the stellar surface: 

dp/p = 0, (21) 



where dp is the Lagrangian pressure perturbation. Not only is this result in agreement with previous results, but it 
also specifies how fast 5p goes to zero near the stellar surface. More details on this method are given in Appendix A. 



4 



D. Reese et al.: Acoustic oscillations of rapidly rotating polytropic stars 



This new choice of variables leads to the following set of equations: 

\b = -Nv ■ VH - HV ■ v, (22) 
\Hv = -H (Vn + V*) + VH (-Nil + "XlHe z x v, (23) 

= A^-^^b. (25) 

If Ti = 7 then N% = and the above system reduces to: 

XNAU = -Nv ■ VH - HV ■ v, (26) 
Xv = - Vn - V* - 2ne z x v, (27) 
= A* - NAH^IL. (28) 

This simplification occurs when the polytropic relation (1) is also the equation of state, a situation typical of white 
dwarfs or neutron stars. Furthermore, both II and b become proportional to the Eulcrian perturbation of the enthalpy, 
thus justifying a posteriori the choice of these variables. As a result, apart from a few multiplicative factors, and the 
lack of a dissipative force, this second set of equations corresponds to those obtained by Yoshida & Eriguchi (1995). 



2.4. Domains and boundary/interface conditions 

In order to complete the eigenvalue problem given by Eqs. (22)-(25), it is necessary to specify a number of boundary 
conditions. The basic requirements are that the solutions remain bounded at the surface and at the centre of the star, 
and that the gravity potential goes to zero at infinity. 

At the centre of the star, the regularity conditions are classically expressed in terms of spherical harmonics (see 
Eqs. (53) and (54)). By using the variables II and b from the generalised Frobenius study, the solution is naturally 
bounded on the star's surface. However, the use of these variables leads to a degeneracy between Eqs. (22), (24) and the 
radial component of Eq. (23) on the surface of the star. This problem is remedied by replacing the radial component 
of Eq. (23) with its radial derivative on the surface. 

It is also necessary to impose a boundary condition on the perturbation to the gravity potential in order to 
ensure that the potential goes to zero at infinity. Traditionally, this is done by doing a harmonic decomposition of W and 
imposing the correct condition on each component. However, such a procedure becomes complicated on a spheroidal 
surface, and it is not certain whether the decomposition of will converge for highly flattened configurations (Hachisu 
et al., 1982). We therefore employ a different method based on Bonazzola et al. (1998). It consists in adding a second 
domain V2 which is bounded on the inside by the star's surface and on the outside by a sphere of radius r = 2 (which 
is twice the equatorial radius). We solve Poisson's equation in this domain and impose the correct boundary condition 
on its outer boundary (where we can safely apply a harmonic decomposition). On the inner boundary, it is necessary 
to use interface conditions which ensure the continuity of and its radial derivative across the stellar surface. 



2.5. Spheroidal geometry 

The next step in the calculations is the choice of a coordinate system based on Bonazzola et al. (1998) for each domain. 
In order to preserve spectral accuracy, the system of coordinates in the first domain needs to fit the surface of the 
star, and provide a non-singular transformation in the centre. As in Paper I and Ricutord et al. (2005), we choose the 
following definition for the radial coordinate £, which ensures a good convergence of the numerical method: 

r(C, 0) = (1 - e)C + ~ (Rs(0) - 1 + e) , (29) 

where e is the flatness given by Eq. (4), (r((, 9), 0, </>) are the spherical coordinates corresponding to the point (C, 0, 4>), 
and R s {0) is the surface of the star. By setting £ = 1, one obtains r(l,6) = R s (0), and the centre r = is given by 

C = o. 

In second domain, we used the following definition: 

r(C, 0) = 2e + (1 - e)C + (2C 3 - 9C 2 + 12£ - 4) (R s (0) - 1 - e) , (30) 

where C € [1, 2] . This mapping is chosen so as to insure the continuity of r and across the boundary £ = 1, and so 
that the surface given by ( = 2 corresponds to the sphere r = 2 (r^ denotes d^r). 
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Figure 1. Coordinate system used in 
computing the equilibrium model of 
the star and the pulsation modes. The 
domain V corresponds to the star itself 
(which in this case is a N = 3 poly- 
trope at 84% of the breakup rotation 
rate). The domain Vi encompasses the 
star, its outer limit being a sphere of 
radius r — 2 (twice the equatorial ra- 
dius). The dotted lines correspond to 
C = 0.1, 0.2. ..0.9 in the domain V, 
£ = 1.1, 1.2 ...1.9 in the domain V2 and 
9 = 0°, 18°, ...342°. The centre of the 
star corresponds to C, = 0, the star's 
surface (the boundary between V and 
V2) to £ = 1, and the outer bound- 
ary to f = 2. The continuity equation, 
Euler's equation, the energy equation 
and Poisson's equation are solved in 
the first domain. The letters v, P 
and p in the domain V show that these 
variables intervene in the first domain. 
In the second domain, only Poisson's 
equation is solved, and only the pertur- 
bation of the gravity potential * inter- 
venes. This is represented by the letter 
$ in the domain V2 ■ 



Once the coordinate system has been established, it is also necessary to choose a set of vectors as a basis. We define 
the following vectors, which are derived from the natural covariant basis (E^, Eg 7 E<p) (defined as Ei = dir): 



a<t> 



c 2 



-Et 



2 



= c! e 

9 c r j 



Eg 



r 1 r( _ 



{rge r + reg) , 



(31) 



\&4> — e <pi 



9 • n V 

r^r^smv rr^ 



where (e r ,eg,e^) are the usual spherical vectors. The vectors (a^,a$, a^) have been chosen so that they become 
(e r , eg, e^) in the spherical limit. Using this base of vectors, we can then express the velocity field as follows: 



v = u^a^ + u e ag + v^a^. 

With these definitions, it is now possible to give an explicit expression of the oscillation equations: 



(32) 



Xb = - 



( 2 Hr c u( (Hrgu 6 

n "T" o 



C 2 r&u c ((r 2 +r 2 )u e 



An- 



A^ 

r C 

Tib 

(JV + 1)A 



r 2 rQ 



Htu C + 



Hgu 9 

c 



C 2 H 
r 2 TQ 



n f 2u c dgu 6 cot Ou 6 dM 



c c 



c 



Csin6> 



X K (r, sjn + r cos^ _ _ + He / b _ \ 
rr^ H \A ) 

2^C 2 sm^ 2n((r e sin 9 + r cos 6)u d^U 



r 



Ar 2 rv \ 7 



H c u c + 



rrQ 

HgU 9 

c 



sin tt sin < 



(33) 
(34) 
(35) 
(36) 
(37) 
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= "-^d^ + c <d< * -J^ fl * + - 2 Ae^ H"~\ (38) 

where: 

cq = -^-3 (2r Ci rerc i e - r 2 r^ - r\rgg + 2rr 2 - rjr^ - r 2 rg cot 6) , (39) 
r r c 

A 90 = d 2 ee + cot 0d e + ~^-d 2 H . (40) 
sin # 

Eq. (33) is the continuity equation, Eqs. (34)-(36) are Euler's equations, Eq. (37) corresponds to the adiabatic 
energy equation and Eq. (38) is Poisson's equation. Euler's equations have been used in their covariant rather than 
contravariant form, as it is advantageous from a numerical point of view. In order to understand this, it is helpful to 
bear in mind that Hg / H converges towards H^g / on the stellar surface, whereas /H is unbounded (since ^ 
on the surface). As a result, the radial component of Euler's equation reduces to b = ANU on the surface (which 
incidentally is already implied by a linear combination of the energy and continuity equations), whereas the two other 
components retain useful information on the surface. If the equations where in their contravariant form, than the 9 
component of Euler's equation would also reduce to b = ANU on the stellar surface since would appear in this 
equation, thus preventing the possibility of dividing by H. It would then be necessary to also replace this equation by 
a supplementary boundary condition which would provide the information already contained in the covariant form. 
This system of equations applies in the first domain, except for Poisson's equation which is used in both domains (of 
course, the density perturbation no longer appears in the second domain). 



2.6. Numerical Method 

In order to solve Eqs. (33)-(38), we project these equations onto the spherical harmonics (Rieutord, 1987). This is 
done in two steps (c.f. Paper I). First of all, the different unknowns are expressed in terms of a sum over the spherical 
harmonics. Explicitly, we obtain: 

oo 

b = E ( 41 ) 

i'=\m\ 
oo 

n = £ i£yF, (42) 

i'=\m\ 
oo 

* = ]T (43) 

i'=\m\ 
oo 

v= J2 uiB% + visy + W iTy, (44) 

t' = \m\ 

where Y™ is the spherical harmonic of degree £' and azimuthal order m and b m , IT^ etc. are radial functions that need 
to be determined, and which only depend on (. The equilibrium model is axisymmetric meaning that the variable 
<f> is not coupled to the two others variables. Therefore, there is no summation over the azimuthal order m in these 
expressions. However, £ and 9 are not separable since the star does not respect spherical symmetry. As a result, it is 
necessary to sum over the harmonic degree £'. 
Rf,Sf, and are defined as follows: 

HP = Y e Ta ( , (45) 
S% = d e Y e Ta e + D^a*, (46) 
T™ = D+YFae - dgY^a^, (47) 

D = (48) 
sin0 

It is worth noting that R™, S™, and T™ are not the usual vectorial spherical harmonics because (a^,ag,a^,) is not 
the same as (e r , eg, e^). However, in the spherical limit, they will become the usual spherical harmonics. An explicit 
expression for each component of the velocity reads: 



^ = Y ™ u i> ( 49 ) 



t' = \m\ 
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e>=\m\ 

oo 
t'=\m\ 



11! , , 



(50) 
(51) 



Once the unknown quantities have been expressed this way, the next step is to project the equations themselves onto 
the spherical harmonic basis. Eqs. (33), (34), (37) and (38) are multiplied by {Yf 1 }* and integrated over Air steradians. 
For each harmonic degree £ of {Yf 1 }* , a different equation is obtained. The remaining equations are obtained from 
Eqs. (35) and (36) in a more complicated manner. We compute the integral over Air radians of {Eq. (35)} {d$Y™}* + 
{Eq. (36)}{D F™}* and {Eq. (35)} {D F™}* - {Eq. (36)} {d g Y e m }* . This operation corresponds to what would be 
done in the spherical case (i.e. a projection onto the vectorial spherical harmonics). As a result, the system thus 
obtained reduces to the classical uncoupled system of equations in the spherical limit. In general, however, this set of 
equations is a highly coupled system of ordinary differential equations in terms of the radial coordinate C, the solution 
of which gives the unknown radial functions (see Appendix B). 

In order to solve this system numerically, we first begin by using a finite number of spherical harmonics L max - 
The equations arc then discrctised onto a Gauss-Lobatto collocation grid of N T + 1 points, based on the Chebyshev 
polynomials. This results in an algebraic system of the form Av = XBv in which A and B are numerically determined 
square matrices. The eigensolutions (A, v) of this system correspond to the frequencies and pulsation modes of the star. 
They are determined iteratively through the Arnoldi-Chebyshev algorithm (e.g. Chatelin, 1988). The coefficients of 
matrices A and B are computed using an equilibrium model with a harmonic resolution L mo d and a Chebyshev (radial) 
resolution of N T + 1. They arc calculated using the coupling integrals given in Appendix B. This is achieved through 
Gauss' quadrature method with L rcs points. Typical values for the different resolutions are: iV r = 60, L mo d = 50, 
Lmax — 80, L rcs = 230. 

At this point, we can write the boundary condition on the gravitational potential and the regularity conditions at 
the centre. The boundary condition is applied along the surface r ext = 2 (or ( = 2) on each harmonic component of 
the gravitational potential perturbation (Hurley et al., 1966): 

1 «~ f + V ro = 0. (52) 



1 - e AC, 1 r ext 

The regularity condition depends on the parity of I in a solution. Thanks to star's equatorial symmetry, modes will 
either be described by a sum of spherical harmonics with even degrees or odd degrees 2 (see Sect. 2.7). For modes with 
even harmonics, we apply the following condition at r = (or £ = 0): 

^=0, ^p=0, ^p=0, «*, = (), t4 = 0, «4 = 0. (53) 
The other modes follow the condition: 

q> e =0 IT"=0 b e =0 ^"=0 ^L = ^"=0 (54) 

d( d( dC 



2.7. Mode classification and symmetries 

A number of useful pieces of information can be deduced from the various symmetries present in the system. These help 
with mode classification, reduce numerical demand and explain certain properties which were observed in perturbative 
calculations. 

The first and most obvious symmetry stems from the fact that the equilibrium model is axisymmetric. This implies 
that modes will have a well defined azimuthal order m (as explained earlier on). A second equally obvious symmetry 
results from the star being symmetric with respect to the equatorial plane. This leads to oscillation modes which are 
either symmetric or antisymmetric with respect to the equatorial plane and which are called even or odd, respectively. 
In terms of spherical harmonics, even modes are made up of harmonic components such that I + m is even, except for 
the toroidal component of the velocity field in which £ + m is odd. Odd modes correspond to the opposite situation. 
From a numerical point of view, eigensolutions are described with half as many components as solutions with no 
particular parity. 

There are two more symmetries which are a little more subtle than the previous ones. The first one only applies 
if the Coriolis force is neglected and only shows up in the rotating frame. In this situation, only even powers of the 

2 The toroidal components w l m have the opposite parity with respect to the other components. 
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rotation rate show up in the equilibrium and pulsation equations. As a result, a given mode will also only depend on 
SI 2 and will be a solution for the rotation rates SI and —St. When this symmetry is combined with the next one, then 
for a given multiplet, modes with azimuthal orders m and — to have the same frequency, as was already pointed out 
in Paper I. 

The last symmetry applies even with the Coriolis force and for both rotating and non- rotating frames. Let us 
consider a solution (u, p, P, v, ^, St, m) (we include the rotation rate and the azimuthal order for the sake of clarity) 
and denote by S the operator which gives the mirror image with respect to the meridian passing through <p = 3 . We 
then find that (uj,Sp,SP,Sv,S^, —SI, —to) is also a solution (this is not to be confused with the previous symmetry 
for which (u, p, P, v, — Si, m) was the corresponding solution). This symmetry was pointed out by Clement (1989), 
however some of its consequences on perturbative calculations were not fully appreciated at the time. Let us consider 
a perturbative description of two frequencies with the same radial order and harmonic degree but with opposite 
azimuthal orders. We will obtain expressions of the following form: 

u n , t , m {Sl) = <^ ro +<^ m + w 2 i£;TO 2 + ... + O(fi fe ), (55) 
u n ,e,-m(Sl) = u J l^_ m + ^l^_ m S} + u;l^_ rn Sl 2 + ... + 0(Sl k ), (56) 

where u) J n e m is the j th perturbative coefficient of oj n ,e,m(Sl)- If we apply the symmetry, we find that 0J n ,e,m(Sl) = 
u> n ,e,-m{-Sl). By equating like powers of w„ i ^ i _ TO (-0) with uj n j tm (Sl), we find that uj 3 n e m = (-1)^^ _ m . 
Therefore u>l „ „ will be an even function of to when j is even and an odd function of to when j is odd. This 
explains why the second order coefficients of Dziembowski & Goode (1992) were polynomials in to 2 , and is also found 
at third order (Goupil, private communication). This symmetry can also be used to increase the accuracy of least 
squares estimates of the coefficients based on non-perturbative calculations (see Sect. 4.1.1). 

3. Analysis of the accuracy of the results 

In order to check whether the results presented here are correct, it is important to do a number of internal tests and 
comparisons with previous studies. We first begin by discussing the accuracy of the underlying polytropic models. This 
is then followed by a series of comparisons with other studies. In the two first comparisons, the previous results have 
a limited accuracy, therefore only allowing a qualitative evaluation. The next two comparisons are with very accurate 
results, thus allowing a quantitative evaluation of the precision of the present results. These are then followed by a 
test based on the variational principle and an analysis of the sensitivity of the results to the parameters used in the 
numerical method. Finally, we conclude by estimating the overall accuracy of the results. 

3.1. Accuracy of the polytropic models 

There are several different tests which give an idea of the accuracy of the polytropic models. One way is by looking 
at the effects of different input parameters, such as the radial or harmonic resolution, on various non-dimensional 
parameters like those in Eq. (4). For non- rotating models, these non-dimensional parameters can be compared with 
those given in Seidov (2004). In Table 1, we give such a comparison, which shows that it is possible to correctly obtain 
6 digits after the decimal point. Table 2 contains a and A for an N = 3 polytrope rotating at 0.59 Six- This table 
shows the strong influence of N T and the need for a sufficient radial resolution. It also suggests a precision of 6 digits 
after the decimal point, if we compare the values for N T = 50 and N T — 60. 

Table 1. Non-dimensional parameters of a non-rotating N = 3 polytrope. I/ max = 50 for all the calculations. It is difficult to 
accurately obtain the 7 th digit after the decimal point, in comparison with the values of Seidov (2004). 



A r 


a 


A 


50 
60 
100 

Seidov (2004) 


54.182 480 87 
54.182 481 06 
54.182 480 87 
54.182 481 11 


47.566 520 74 
47.566 520 85 
47.566 520 74 
47.566 520 88 





3 In spherical coordinates, S is defined as follows for a scalar quantity: SA(r, 9, <j>) — A(r, 9, —<j>). For a vector field it takes on 
the definition: SV(r, 9, <jj) — V r (r, 9, —(f>)e r + Vg(r, 8, —4>)e$ — V^(r, 8, — 0)e^. 
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Table 2. Same as Table 1 but for an N = 3 polytrope rotating at 0.59 fix- The radial resolution N T has a stronger effect on 
the values of a and A than the harmonic resolution Z/ max . 



N T imax « A 



50 


16 


81.108 265 69 


63.025 583 86 


20 


50 


81.108 444 82 


63.025 591 55 


50 


50 


81.108 249 13 


63.025 575 39 


50 


60 


81.108 249 08 


63.025 575 36 


60 


50 


81.108 249 38 


63.025 575 52 


60 


60 


81.108 249 38 


63.025 575 52 



In addition to the previous test, it is also possible to apply the virial theorem to obtain a measure of the accuracy 
of the model's structure. In what follows we use the following formulation of the theorem: 

o = J Po ny sin 2 edv+ l -j^ Po ^ dv + J p a dv, (57) 

where p = H N , P Q = H N+1 and $f = tp /h c . For a sufficient number of iterations, it is possible to attain a precision 
of 10~ 13 on the virial test. Beyond this point, successive iterations are useless and can actually decrease the accuracy 
of the model. In Fig. 2, we follow the evolution of A and the virial error with each iteration. As can be seen in the 
figure, there are two phases: a first phase in which the model is approaching the mathematical solution to the problem, 
and a second phase in which the maximum precision has been attained and the model is slowly drifting towards less 
accurate solutions. For some of the rotating configurations and with a well adjusted resolution, this second phase does 
not contain a slow drift but remains close to a fixed point. Either way, the best point at which to stop the iterative 
scheme is at the transition between the two phases. 



)xl0" 7 r Seidov (2004) 



8x10" 



o 

o 7xl0- 7 
to 

CO 



it 6xl0- 7 ~- 
I 



5x10" 



4x10" 



Nr=60 



Nr=50 




-4X10" 11 



-2X10" 11 
Virial test 



Figure 2. A plot of the value of A as a function of 
the virial error. Each iteration is marked by a plus 
"+" and connected consecutively. As shown in the 
figure, the iterated models reach a point of closest 
approach to the mathematical solution, and then 
slowly drift towards less accurate models. 



The models on which are based the pulsation frequencies do not attain as high a precision, because the iterative 
scheme was stopped before the transition between the two phases. This is because we use a small parameter called e 
which controls the relative error on the enthalpy and serves as a stopping criteria. If the value of e is too low, than 
the iterative program never reaches this precision on the enthalpy and therefore does not output the stellar model. 
We therefore set e = 1CT 8 in most calculations, which ensures successful convergence but reduces the accuracy of the 
model. As a result, the virial test typically attains a precision of 4 x 10~ 10 . For the non-rotating model, a takes on a 
value around 54.182473, which starts differing at the 5 th digit after the decimal point from the value given in Seidov 
(2004) and corresponds to a relative precision of ~ 10~ 7 . 



3.2. Comparison with Saio (1981) 

Saio (1981) gives second order perturbative calculations for polytropic models. Based on his coefficients, it is possible 
to obtain pulsation frequencies via the following formula: 

w = wo - (1 - Ci) mCl +{(X 1 +X 2 + Z) + m 2 (Y 1 + Y 2 )} — . (58) 
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In order to compare our results with his, it is necessary to extract perturbative coefficients from our own results. The 
procedure used to find these coefficients is fully described in Sect. 4.1.1. Before applying this procedure, we first had 
to express our results in the same units as Saio (1981) via the following conversion rule: 

w S 8i = \/3aww, (59) 

= \f^^ (60) 

where the subscript "nr" means "non-rotating" (i.e. the value of the parameter for the non- rotating polytrope) and the 
subscript "S81" means that the quantity is in Saio's units. It turns out that in Saio's units, the mass of the polytrope 
depends on the rotation rate. A comparison between his coefficients (u; 2 , Ci, X\ + X2 + Z,Y\ + Y2) and ours showed 
a qualitative agreement between the two (to within 2 %). This reduces the possibility of programming errors affecting 
our results. 



3.3. Comparison with Clement (1984) 

Further comparisons can be done with Clement (1984), who applied non-perturbative techniques to calculate pulsation 
frequencies for N = 1 and N — 3 polytropes. His frequencies are given in the same units as ours and no conversion is 
needed. However, Clement (1984) used a rotational parameter which he called "a" (denoted here as acs4 so as to avoid 
confusion with a from Eq. (4)) and which is "neither dimcnsionless nor scale-free". Therefore, based on the conversion 
given in Tables 1 and 2 of Clement (1984) and following his recommendations, we used the parameter v — Q 2 /2irGp c 
instead. We allowed for uncertainties in the last digit off and therefore calculated a corresponding range of frequencies. 
For example, if v = 1.69 x 10~ 2 we would calculate the frequencies corresponding to v = 1.69 x 10~ 2 ±5 x 10~ 5 . The 
results are presented in Fig. 3. 



1.5 



1.0 
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0.00 
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Clement (1984) 



0.02 0.04 

v = Q 2 /2uGp c 
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Present work 
Clement (1984) 



0.001 0.002 

v = Q 2 /27TGp c 



0.003 



Figure 3. A plot showing the frequencies in Clement (1984) and our calculations of the same frequencies 



As can be seen from Fig. 3 the two sets of results qualitatively agree, which once more makes numerical programming 
mistakes unlikely in our calculations. However, Clement's results usually do not lie in the frequency intervals we calcu- 
lated (this would require a 4-digit accuracy). This is partially due to the fact that it is difficult to accurately reproduce 
the polytropic models he used. In order to illustrate this, we can use the different parameters (v, p c j p, R cq / R po \, 9e/ 9p) 
provided in Tables 1 and 2 of Clement (1984), allow for uncertainties in the last digit, and calculate the corresponding 
ranges for ft*. For a given rotation rate, if all the digits in the four parameters are accurate, then the four different 
ranges for ft* should overlap and give a more precise idea as to the underlying model. However, it was only possible to 
obtain at most three overlapping ranges, and not four. A typical example for the N = 1 polytrope (with ac84 = 0.004) 
is: 



v = 3.57 x 10~ 2 ± 5 x 10~ 5 => ft* = 0.4615 ± 4 x 10~ 4 , 

p c /p =3.40±5xl0~ 3 ft* = 0.4580 ±9.4 x 10~ 3 , 

R eq /R po \ = 1.162±5xl0- 4 ft* = 0.4590 ±7 x 10" 4 , 

9e /g p = 0.686 ±5 x 10~ 4 => ft* = 0.4621 ±4 x 10~ 4 . 



(61) 
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This shows that the error bars we used on Clement's parameters are too small and that the uncertainties on his models 
are larger. Nonetheless, these uncertainties do not fully account for the discrepancies between our frequencies and his. 
This can be shown by the fact that even for non-rotating configurations (where there is no ambiguity on the underlying 
model) the differences are still of the same order of magnitude. 

In general, a quantitative comparison between our results and those of Saio (1981) and Clement (1984) showed 
an agreement to 2 or 3 significant digits. While providing a correct qualitative picture, the precision of these studies 
is insufficient for future missions such as COROT. COROT will observe pulsation frequencies within the range of 
0.1 — lOmHz with an accuracy of 0.6 /xHz for the 20 day runs and 0.08 pHz for the 150 day runs (e.g. Baglin et al., 
2002), meaning that an accuracy of 3 to 5 digits is required. Therefore, it is important to show that our results meet 
up to this requirement through other more constraining tests and comparisons. 

3.4. Comparison with Christensen-Dalsgaard & Mullan (1994) 

Christensen-Dalsgaard & Mullan (1994) give very accurate frequencies for several non-rotating polytropic models. 
Comparing our results with theirs provides a robust test for accuracy. In order to convert our frequencies ui into their 
units, we apply the following conversion rule: 

vcdm = VgViau, (62) 

where j/cdm is our frequency in their units, v g = 99.855377 /xHz, and a is given by Eq. (4). A comparison between 
their frequencies and ours revealed a very good agreement (Auj/uj <~ 10~ 7 for a N = 3 polytrope and Auj/uj <~ 10~ 8 
for a N = 1.5 polytrope at CI = 0). The modes which were compared are: I = to 3, n = 1 to 10 for A" = 3 and 
n = 15 to 25 for N = 1.5. The differences come from round-off errors in the last digit (if we keep the same number of 
significant digits). 

3.5. Comparison with Paper I 

We finally compared our results with those of Paper I. In order to do this comparison, it is necessary to remove the 
Coriolis force and to make the Cowling approximation. No conversion rule is necessary since both sets of results are 
given in the same units. The two sets of frequencies agree quite well, even at large rotation rates (Auj/uj ~ 10~ 7 ). This 
result is significant due to the fact that the set of equations used in Paper I is entirely different than the one used here. 

3.6. Variational test 

The variational principle provides an integral formula which relates a pulsation frequency to the structure of the 
corresponding mode (e.g. Lynden-Bell & Ostriker, 1967). It is therefore possible to apply this formula to a numerically 
calculated eigenmode to obtain a "variational frequency". The error on this frequency is quadratic in the error of 
the eigenfunction (c.f. Christensen-Dalsgaard & Mullan, 1994). By comparing this frequency with the one obtained 
directly, it is possible to estimate the accuracy of the results. We used the following non-dimensional formulation of 
the variational principle: 

u 2 f Po \v\ 2 dV - M 2 ( f - f I V*| 2 dv) + 2iu> f Po n ■ (v* xv)dV- f Po N 2 \v ■ e g \ 2 dV = 0, (63) 

Jv \JV PoC Q J Voo J Jv Jv 

where V is the volume of the star, V^> is infinite space, and e g the unit vector in the same direction as the gravity 
vector. The pressure is replaced by H N H and the different integrals are performed numerically 4 which gives a second 
degree equation in uj. Solving this equation gives the variational frequency which can then be compared with the direct 
calculation of ui. Generally, we find differences Auj/uj ~ 10~ 8 or better between the two. This can be compared with 
the results of Ipser & Lindblom (1990) who found differences of 10 -3 when they applied the variational principle to 
their calculations. An explicit formulation of the variational principle in spheroidal geometry is given in Appendix C. 

3.7. Influence of the parameters from the numerical method 

A final test consists in modifying different input parameters and seeing the effect it has on the results. We have 
therefore applied this test to a few modes which are representative of all the modes that have been calculated. The 
parameters that were modified are the radial resolution N T (which is the same for the equilibrium model and the 



4 The numerical integration was based on Gauss' quadrature method and a spectral expansion, using a radial resolution of 
101 points and an angular resolution of 200 points. 
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pulsation mode), the harmonic resolution of the model L mo d, the harmonic resolution of the pulsation mode I/ max , 
the shift a used in the Arnoldi-Chcbyshcv algorithm 5 , and e which controls the relative error of the enthalpy in the 
equilibrium model. Table 3 lists the values used for the different parameters and the induced frequency variations. 
For a given parameter, we used the frequency obtained at highest resolution (or lowest value for e) as a reference. In 
most cases, we obtained a rough plateau at the different levels given in Table 3. In some cases however, there was a 
definite decrease of the error. For instance, for those modes in which it was tested, the error was roughly proportional 
to e. Also, for high frequency modes, the error strongly decreased as N r increased, as could be expected for high radial 
orders. In the table, we put the lower/final values of the error for both of these parameters. The information on the 
shift is slightly different. The line "Values" gives the amplitude of the variation on the value of the shift. The next 
two lines contain the standard deviation of the results. 

The results on e are not representative of the calculated frequencies. As was pointed out in Sect. 3.1, the number 
iterations was usually less than optimal because of a large value of e (1CT 8 instead of 1CT 10 ), thus resulting in a 
decreased accuracy. The relative error on low frequency modes is 1CT 8 and that of high frequency modes IO -7 . 

Table 3. Frequencies variations in terms of different parameters. L max is the harmonic resolution of the pulsation modes, L mot j 
the harmonic resolution of the models and N T the radial resolution. The line "Values" gives the different values that were used 
for the resolutions, e and the width giving the variation of the shift. The two following lines give the order of magnitude of the 
induced frequency variations (in units of ^4mGp c ) ■ 







^mod 


N T 


Shift 


e 


Values 


40, 44 ... 80 


30, 34 .. 


. 70 32, 36 .. 


.60 2 - 


5 x 10~ 4 


io- 8 ...io- IU 


Low frequency modes 


< IO -15 


io- 10 


1Q -io 


10" 


13 


10 -io 


High frequency modes 


1Q -14 


io- 9 


10" 10 


10" 


■11 


10" 9 



3.8. Discussion 

Overall, the main source of error in the present calculations is the uncertainties on the equilibrium model. This is 
because we chose a convergence criteria which was sure to be met, but which lead to a number of iterations less than 
optimal. This therefore leads to a global accuracy of 7 digits after the decimal point (in units of \/4irGp c ), the last 
digit being uncertain. Table 3 however shows that these calculations could potentially be made more accurate. The 
present accuracy is nonetheless largely sufficient for the requirements of COROT, which will be at most 5 significant 
digits. 

4. Results 

We now proceed to present the results themselves. We followed acoustic adiabatic pulsation modes (with Ti — 5/3) 
from a zero rotation rate to 0.59Qk (where Qk — \J GM / is the Keplerian break-up rotation rate), using the same 
procedure as in Paper I. This involves identifying the frequencies at Q — 0, following their evolution while progressively 
increasing the rotation rate, and working through a number of avoided crossings. The underlying polytropic models 
have an index N = 3 which gives a polytropic exponent 7 = 4/3. The modes that were calculated are: £ = to 3, 
n = 1 to 10 and m = — £ to I both with and without the Coriolis force. 

4.1. Comparison with perturbative methods 

In this section, we compare complete and perturbative calculations so as to determine the range of validity of pertur- 
bative methods. 

4.1.1. Perturbative coefficients 

In order to compare perturbative calculations with complete ones, it proved necessary to compute our own perturbative 
coefficients, since we were unable to find perturbative coefficients for polytropic models with a sufficient accuracy in the 
literature. Instead of using the traditional method of perturbing the fluid equations and finding corrections of various 

5 The shift comes from shift-and-invert methods and corresponds to a trial value around which the Arnoldi-Chebyshev 
algorithm looks for frequencies. See Valdettaro et al. (2006) for an extensive discussion on the role of the shift in numerical 
errors. 
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orders on the frequencies (see Soufi et al. (1998) for a complete description), we did a series of complete calculations 
for small rotation rates = 0, 1CT 6 , 1CT 5 , 1CT 4 , 1CT 3 , 0.002, 0.004 ... 0.018) and applied a least squares fit to 

the results. In order to increase the accuracy of such calculations, we made use of Eq. (55) and separated even and 
odd powers of the rotation rate: 

T '" +< {>m fi 2 + < {>ro fi 4 + e>(fi 6 ), (64) 



2 



UnXm 2 Un>l -- m - < l , fB n + < Iira n» + < Iim n 6 + o(n 7 ). (65) 

By fitting {u) n ,l,m+<^n,l, -m)/2 and (u) n j, m —u n ,i, -m)/2 the number of unknowns is reduced to three and the residues 
are smaller. It is necessary to include the fourth and fifth powers of the rotation rate so as to ensure that the second 
and third order coefficients are reasonably accurate. The results are given in Table 4 for frequencies and rotation rates 
in units of From these coefficients the frequencies are given through the following formula: 

w = u} -m(l- C)Q + (Di + m 2 D 2 ) O 2 + m {T x + m 2 T 2 ) n 3 + O (ft 4 ) . (66) 

The form of the second degree coefficients was obtained from Saio (1981) and that of the third degree coefficients from 
Goupil (private communication). In order to express these results in units of Qk instead, one can use the following 
perturbative formula: 

Al^LWoll") >, (67, 



where tt p / = (GM/R^ and A ~ 0.77166. 



Table 4. Perturbative coefficients for a N = 3 polytrope, deduced from complete calculations. The frequencies and the rotation 
rate are expressed in units of (GM/H^ oi ^ 1 ^ 2 . 



wo 


C 


Di 


D 2 


Ti 


LOO 


C 


Di 


D 2 


Tj 


T 2 






£ = 










1 = 2 








3.042155 




-1.194 






3.906874 


0.1538359 


-1.578 


-0.294 


0.02344 


0.00527 


4.121230 




-2.315 






5.169469 


0.0818188 


-2.396 


-0.459 


0.00493 


0.00477 


5.336900 




-3.439 






6.439990 


0.0544285 


-3.146 


-0.606 


0.00020 


0.00307 


6.591212 




-4.484 






7.708951 


0.0403695 


-3.867 


-0.745 


-0.00087 


0.00218 


7.855027 




-5.484 






8.975891 


0.0318248 


-4.572 


-0.880 


-0.00094 


0.00169 


9.120432 




-6.459 






10.240946 


0.0260651 


-5.267 


-1.013 


-0.00074 


0.00139 


10.384948 




-7.416 






11.504260 


0.0219090 


-5.955 


-1.144 


-0.00049 


0.00119 


11.647767 




-8.361 






12.765953 


0.0187647 


-6.636 


-1.273 


-0.00025 


0.00105 


12.908679 




-9.298 






14.026134 


0.0163027 


-7.314 


-1.402 


-0.00006 


0.00094 


14.167704 




-10.228 






15.284901 


0.0143246 


-7.987 


-1.530 


0.00010 


0.00086 






l = \ 










£ = 3 








3.377036 


0.0295367 


-1.072 


-1.030 


-0.04612 


4.294602 


0.1193654 


-1.898 


-0.169 


-0.01155 


0.00041 


4.642432 


0.0342809 


-1.760 


-1.699 


-0.04204 


5.591067 


0.0742468 


-2.728 


-0.240 


-0.00113 


0.00157 


5.909240 


0.0335303 


-2.402 


-2.315 


-0.02715 


6.878680 


0.0517251 


-3.496 


-0.307 


0.00015 


0.00140 


7.176668 


0.0305143 


-3.019 


-2.901 


-0.01634 


8.158826 


0.0387755 


-4.236 


-0.371 


0.00038 


0.00113 


8.443277 


0.0270732 


-3.621 


-3.470 


-0.00951 


9.433911 


0.0305232 


-4.960 


-0.434 


0.00043 


0.00091 


9.708372 


0.0238467 


-4.213 


-4.028 


-0.00524 


10.705348 


0.0248710 


-5.674 


-0.496 


0.00044 


0.00076 


10.971700 


0.0210064 


-4.797 


-4.578 


-0.00254 


11.973956 


0.0207887 


-6.380 


-0.557 


0.00046 


0.00064 


12.233222 


0.0185621 


-5.375 


-5.122 


-0.00080 


13.240238 


0.0177183 


-7.081 


-0.618 


0.00048 


0.00055 


13.492998 


0.0164731 


-5.950 


-5.662 


0.00034 


14.504529 


0.0153345 


-7.777 


-0.678 


0.00050 


0.00048 


14.751133 


0.0146882 


-6.521 


-6.198 


0.00108 


15.767068 


0.0134359 


-8.470 


-0.738 


0.00051 


0.00042 



In order to estimate the accuracy of these perturbative coefficients, there are a number of tests that can be done. 
First of all, the zeroth order coefficients are simply the pulsation frequencies without rotation but are treated as 
unknowns in the least squares development. The frequencies without rotation are recovered in the least squares fit to 
an accuracy of at least 5.4 x 10~ 8 in the units of Tabic 4. The first order coefficient C, can be calculated via integrals 
based on the zeroth order solution (Ledoux, 1951). This alternate way of calculating the coefficients agrees to within 
1.4 x 10~ 9 (this does not necessarily mean that the coefficients are accurate to that precision but does show a high 
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degree of internal coherence). For the second and third degree coefficients, we checked to see if they satisfied the 
forms given in Eq. (66); the number of significant digits in Table 4 has been adjusted accordingly. These forms were 
a constraint only on the I = 2 and 3 second order coefficients and on the I = 3 third order coefficients. Another test 
we did consisted in applying the least squares fit to a subset of the results used in the first fit and seeing whether the 
coefficients were altered. This test indicates roughly the same accuracy as the other tests. 

4.1.2. Comparison 

Based on these coefficients, it is possible to calculate perturbative frequencies which can then be compared to the 
complete calculations, thereby establishing a domain of validity for perturbative methods. In Fig. 4, we show two 
such domains for 3 rd order methods, one for each of COROT's error bars (0.6 /iHz for the 20 runs and 0.08 /iHz 
for the 150 day runs). The underlying polytropic models have a fixed mass of 1.9 Mq and a fixed polar radius of 
2.3 Rq, both of which are typical of 5 Scuti pulsators. When the distance between the perturbative frequency and the 
complete one exceeds COROT's error bars, the frequency is shown in black. Otherwise, it is shown in grey. From these 
figures, it is clear that complete methods are required beyond v ■ sini = 75km.s _1 for COROT's 20 day programs and 
= 50km.s _1 for COROT's 150 day programs. 
It is important to bear in mind that the domain of validity obtained for perturbative methods depends on the choice 
of rotational variable used in the development. In order to illustrate this, suppose we develop a frequency in terms of 
two different rotational parameters X and Y: u> = a + a x X + a 2 X 2 + a 3 X 3 + 0(X 4 ) = b a + biY + a 2 Y 2 + a 3 Y 3 + 0(Y 4 ). 
When the relationship between X and Y is more complex than a simple proportionality, the neglected terms, 0(X i ) 
and 0(Y 4: ), are not the same. As a result, a 3 rd order development in terms of X or Y will give different values for 
u>, thus modifying the corresponding domain of validity. Therefore, we decided to compute the domain of validity 
associated with the variables Q/Qk and Q/ (GM / R 3 ol ) 1/>2 to see if there was a substantial difference between the two. 
For individual frequencies, there can be large differences, but when all the frequencies are considered, the global result 
is roughly the same. 

In Fig. 5, we show the differences between complete frequencies and perturbative ones at 0.59£1k- We have kept 
the same parameters for the equilibrium model as in Fig. 4. As can be seen from the figure, differences between the 
two sets of calculations are substantial and comparable to the large frequency separation (which seems to survive 
rotation). The order of frequencies is not the same between the two sets of calculations. As a result, it is necessary to 
use complete calculations in order to correctly interpret a pulsating star rotating at such a high rotation rate. 

Recently, Suarez ct al. (2005) attempted to model Altair through asteroseismology. The effects of rotation were 
included in the pulsation modes using 2 nd order perturbative methods. Later interferometric studies suggested an 
equatorial velocity of 280 km.s -1 (Domiciano de Souza et al., 2005), which is above 216 km.s -1 , the equatorial velocity 
corresponding to Fig. 5 (if we use a mass of M = 1.8 Mq and a polar radius i? po i = 1.7 Rq instead, we obtain 
v cq = 244km.s _1 ). As a result, it is pretty obvious that what is required in Suarez et al. (2005) is complete calculations 
of the effects of rotation before being able to interpret Affair's oscillation spectrum (not to mention complete models 
of rapidly rotating stars) . 

4.1.3. Relative importance of the Coriolis and centrifugal forces 

It is then interesting to analyse what is the main source of differences between perturbative calculations and complete 
ones. In Fig. 6 we show three different graphs which give the relative errors associated with different calculations: 



where the subscripts (a), (b) and (c) correspond to the different panels in Fig. 6, the subscript "pert." to 3 rd order 
perturbative calculations and the superscript "no Cor." to calculations done without the Coriolis force. From these 
panels, it is possible to deduce the dominant role of the centrifugal force in the differences between perturbative and 
complete calculations. Panels (a) and (b) are very similar, yet the first one includes the Coriolis force and the second 
one excludes it. Panel (c) shows the errors which come from excluding the Coriolis force. These errors are at least ten 
times smaller than in cases (a) and (b) and decrease with the radial order. This decrease is expected because as the 
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Figure 4. Plots of the evolution of pul- 
sation cyclic frequencies {v = ui/2n) 
as a function of the rotation rate. The 
frequencies are followed from the non- 
rotating case to 0.59 Q,k using the pro- 
cedure described in Paper I. At small 
rotation rates, it is easy to recognise the 
usual multiplet structure as predicted 
from perturbative methods. Once the 
rotation rate is sufficient large, the mul- 
tiplets are less regular and overlap, 
which greatly complicates the inter- 
pretation of the oscillation spectrum. 
Superimposed is the domain of validity 
of 3 rd order perturbative methods us- 
ing COROT's error bars of 0.6 /iHz (up- 
per panel) and 0.08 ^iHz (lower panel). 
The calculations were done with N = 3 
polytropic models with i? po i = 2.3-R© 
and M = 1.9M Q . 
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radial order n of the mode increases, the time scale of the oscillations decreases and becomes much shorter than the 
1/fi time scale associated with the Coriolis force. As a result high order modes are less affected by the Coriolis force. 

The effects of the centrifugal force, on the other hand, increase with radial order. The reason for this, as explained 
in Paper I, is that changes in the stellar structure and the sound velocity profile causes modifications which are roughly 
proportional to the frequencies. For spherically symmetric changes in a star's structure, we have A In lo = —A In dr/c 
(where A means the variations due to the change in stellar structure), based on Tassoul's asymptotic formula (Tassoul, 
1980). The same principle applies for more complicated changes in the structure, such as those provoked by the 
centrifugal force, but will have a more complicated mathematical formulation. One way of illustrating this is by 
plotting the ratios D\/u)o (see Eq. (66)), which correspond to dlnoj/dfl 2 calculated at f2 = for m — modes. We 
take the derivative with respect to tt 2 since the effects of the centrifugal force begin at 2 nd order in fi. In Fig. 7, we 
can see that these ratios approach constant values as n increases for a given I. Furthermore, we have plotted these 
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Figure 5. This figure shows a compar- 
ison between perturbative cyclic fre- 
quencies and complete ones at 0.59£Ik- 
The equilibrium model is the same as 
in Fig. 4. The complete calculations 
are represented by the bars that go to 
the right and the perturbative ones by 
the bars going to the left. The length 
of the bars gives the azimuthal or- 
der. Complete and perturbative calcu- 
lations are connected by dotted lines 
(the correspondence is based on mode 
labelling as described in Paper I). It 
is clear from this figure that pertur- 
bative calculations lead to substan- 
tial error at high rotation rates, and 
cannot correctly anticipate the order 
of the modes. Even if the perturba- 
tive frequencies were multiplied by a 
global corrective factor, the agreement 
remains poor. 

Harmonic degree I 

ratios both with and without the Coriolis force, thereby demonstrating that this effect is entirely due to the centrifugal 
force. For non-axisymmetric modes, the relevant ratios {D\ + m 2 D 2 )/wo also show the same behaviour. This shows 
that the effects of the centrifugal distortion is roughly proportional to the frequency. 

It can then be expected that making errors on the effects of the centrifugal force will also lead to differences 
proportional to the frequencies. If we look at panel (b) of Fig. 6 in which the Coriolis force has been suppressed, we 
can see that the relative differences between perturbative and complete calculations Slo jio actually increase with radial 
order (at least for lower and moderate rotation rates). It is an open question whether or not these relative differences 
will approach an asymptotic limit as the radial order increases like in Fig. 7. 

The increase of perturbative errors with radial order may have important implications for stars which pulsate in 
high radial overtones. Examples of these are solar type stars which typically pulsate with radial orders between 15 
and 25. If we assume that perturbative errors scale with frequency and frequency with radial order, we can estimate 
at what point complete methods will be necessary to calculate the effects of rotation on such modes. The average limit 
for n — 25 pulsation modes in M = 1M Q , i? po i = 1 Rq stars would be v cq = 25km.s~ 1 (45km.s _1 ) for COROT's 
primary (secondary, resp.) program. This implies that non-perturbative effects of rotation could be visible for moderate 
rotation rates. Nonetheless, direct comparisons between perturbative and complete calculations for high order modes 
are necessary to confirm this conclusion. 

In order to further understand the role of the centrifugal force in perturbative errors, it is helpful to bear in mind 
the approximations which result from applying perturbative methods. First of all, the frequencies, the stellar structure 
and the mode structure are all described by low degree polynomials in Q. This then leads to the following effects: the 
equilibrium structure is only described by £ = and I = 2 spherical harmonics (for 2 nd and 3 rd order methods) and 
the pulsation mode structure is also limited to a few spherical harmonics. In order to analyse the effects of using only 
a few spherical, we did highly truncated numerical calculations, in which the equilibrium model has been reduced to 
2 spherical harmonics and the pulsation modes to 4 spherical harmonics (2 poloidal + 2 toroidal). In Fig. 8, the black 
lines represent the relative differences between these truncated calculations and complete ones. These differences are 
significant, thus showing the need for more spherical harmonics. Nonetheless, we also plot in Fig. 8 the perturbative 
errors which are higher than the truncated calculations. This shows that including higher order terms in f2 in the 
contribution of even the lowest degree spherical harmonics can improve results. 

4.2. Discussion 

As can be seen from previous sections, differences between perturbative calculations and complete ones can be quite 
substantial. This is problematic because obtaining accurate results is crucial in asteroseismology. Differences between 
theoretical calculations and observed frequencies need to come from differences between the stellar model and the 
star's actual structure rather than from an approximate treatment of the effects of rotation. Otherwise, modifying 
the stellar structure so as to match the observed frequencies will end up compensating errors in the calculation 
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Figure 6. Relative error of different calculations of I = 2, 
m = — 1 modes of various radial orders. The radial orders 
are indicated on the right side of each panel. Panel (a): 
the relative error from using 3 rd order perturbative meth- 
ods. Panel (b): Same as panel (a) except that the Coriolis 
force has been excluded from both perturbative and com- 
plete calculations. Panel (c): the relative error which comes 
from neglecting the Coriolis force. Explicit expressions for 
the various errors are given in Eqs. (68)-(70). The similarity 
between panels (a) and (b) and the differences with panel 
(c) show the dominant role of the centrifugal force in errors 
related to perturbative methods. 



of frequencies instead of improving the stellar model. Moreover, large errors on frequency calculations can lead to 
erroneous mode identifications, especially if proximity between observed and theoretical frequencies is the only criteria 
used for establishing such identifications. Interestingly, establishing a correct mode identification is one of the key 
difficulties in interpreting S Scuti stars (e.g. Goupil et al., 2005). Mode misidentification can occur because frequencies 
are not in the same order in perturbative calculations as in complete ones (in Fig. 5, this can be seen by the dotted lines 
for the I = 1, 2 or 3 modes which cross each other, thus indicating an exchange of position between two frequencies) and 
because one can no longer rely on the usual frequency patterns used in slowly rotating stars (see Paper I) . An erroneous 
mode identification then leads a false interpretation of pulsation frequencies due to an incorrect understanding of the 
geometry of the pulsation mode and the stellar regions which it probes. This problem is further aggravated by the 
fact that perturbative methods only give an approximate idea of the structure of a given mode anyway. Fully taking 
into account the effects of rotation on stellar pulsation increases the likelihood of obtaining a correct identification 
and gives a better understanding of mode structure, especially when the rotation rate is high. However, in order to 
obtain such a mode identification, the underlying stellar model needs to be sufficiently close to reality so as enable a 
successful matching between theoretical predictions and observations. It is possible that even then, mode identification 
is uncertain due to multiple solutions which fit a set of observed frequencies. 



5. Conclusion 

In this work we have explored some of the effects of rapid rotation on stellar acoustic pulsations. This was achieved 
thanks to a numerical method which combines spheroidal geometry and spectral methods, as in Paper I. Through a 
detailed analysis, we have shown that our results have a 6 to 7 digit accuracy. This analysis included a discussion on the 
accuracy of the underlying polytropic models, comparisons with Saio (1981), Clement (1984), Christensen-Dalsgaard 
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Figure 7. Ratio of 2 nd order perturbative coefficients to cor- 
responding frequencies at Q = 0, as a function of the radial 
order n. As n increases, these ratios approach a constant 
value, which shows the proportional effects of the centrifu- 
gal deformation of the star. 



Figure 8. A comparison between relative errors in pertur- 
bative calculations and those in truncated calculations (in 
which the equilibrium model is reduced to 2 spherical har- 
monics and the pulsation mode to 2 poloidal + 2 toroidal 
harmonics). This figure shows that truncated calculations 
are closer to complete ones than perturbative calculations, 
which shows the importance of retaining higher powers of fi 
in the lower spherical harmonics describing the stellar struc- 
ture and mode structure. The uneven aspect of the lines 
corresponding to the truncated calculations are due to in- 
terpolation errors. 



& Mullan (1994) and Paper I, a test based on the variational principle and some tests on the sensitivity of the results 
to the parameters of the numerical method. 

In the future, the satellite COROT is expected to measure stellar oscillation frequencies with a precision of 0.08 ^Hz 
(primary targets) or 0.6 /zHz (secondary targets). In the frequency range considered in this paper, we find that for a 
M = 1.9 Mq, R = 2.3 Rq star, perturbative methods cease to be valid for COROT's primary (secondary) program 
beyond v ■ sini = 50km.s _1 [v ■ sini = 75km.s _1 , resp.). At a rotation rate of 0.59fix, the perturbative spectrum is 
very different from the one based on complete calculations. Therefore, any attempt to interpret stellar pulsations using 
the perturbative approach at comparable rotation rates is likely to fail. Using complete methods on the other hand 
increases the likelihood of obtaining a correct mode identification, and gives an accurate description of the structure 
of pulsation modes. Both of these are crucial when interpreting observed pulsation modes. 

Further investigation has shown the dominant role of the centrifugal force in modifying the frequency spectrum and 
causing perturbative errors. This is because while the effect of the Coriolis force decreases as the frequency increases, 
the effect of the stellar deformation increases roughly proportionally to the frequencies. Therefore, the errors which 
arise from perturbative descriptions of the centrifugal distortion are also amplified in higher order modes. As a result, 
it may be necessary to use complete methods for moderately rotating stars which exhibit high order modes. 

Some of the issues which were discussed in Paper I and have yet to be discussed for the present results include: an 
analysis of the regularities in the oscillation spectrum at high rotation rates and a study of the visibility of the different 
modes based on their structure. These will be the subject of forthcoming papers. A few preliminary examinations have 
already confirmed some of the conclusions given in Paper I, such as a strong equatorial concentration of mode structure 
at high rotation rates, or the transition from one frequency spectrum organisation to another. 

Future work includes working with more realistic models, and studying gravity modes in spheroidal geometry. The 
transition to more realistic models is essential before being able to compare theoretical frequencies with observations. 
Coming up with realistic models that fully include the effects of rotation and in particular the centrifugal distortion 
is no easy task, but is the subject of active research {e.g. Roxburgh, 2004, Jackson et al., 2005, Rieutord et al., 2005). 
Calculating the associated pulsation frequencies and comparing them to observations will provide crucial information 
on stellar structure and enable a better adjustment of these models. 

The study of the effects of rapid rotation on g- modes is of interest for the interpretation of 7 Doradus stars, which 
are g-mode pulsators and can be rapid rotators. Previous studies on the non-perturbative effects of the Coriolis force 
on g- modes (Dintrans et al., 1999, Dintrans & Rieutord, 2000) have revealed their important role in altering the 
geometry and frequencies of these modes. This behaviour is entirely different from that of the high frequency acoustic 
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modes presented here. It is then interesting to understand what the effects of the centrifugal force will be on g-modcs 
and how it will compare with the effects of the Coriolis force. 
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Appendix A: "Generalised" Frobenius method 

A.l. Description 

The starting point in this method is the following equation: 

dY{ Y ,Z) + -M*, V, z)Y(x, y, z) = 0. (A.l) 
dx x 

This equation looks very much like a first order Frobenius equation except that two other variables, y and z, intervene 
in the different terms (for a description of the more traditional version of the Frobenius method, see Bender & Orszag, 
1978). The quantity Y(x, y, z) can be a scalar or a vector. The operator A(x, y, z) can include derivatives in the y and 
z directions and needs to be analytic in the x direction, so that we can write: 

oo 

A(x,y,z) = J2My,z)x n - (A.2) 

n=0 

We then look for the behaviour of Y(x, y, z) along the boundary x = 0. If we develop Y(x, y, z) in the following 
manner, 

oo 

Y(x,y,z)-]TY„( y ,z)x"+ c \ (A.3) 

n=0 

then we obtain the following zeroth order equation: 

aY (y,z) + A (y,z)Y (y,z) = 0, (A.4) 

where a is the leading power of x in Y(x,y, z). Therefore, to obtain a, one needs to solve an eigenvalue problem in 
terms of the coordinates y and z, along the entire surface x = 0. The remaining Y n are defined through the following 
recurrence relation: 

71-1 

n>l, [{a + n)Id + A ]Y n = -J2 A n-kY k , (A.5) 

where Id(Y) = Y. This series is defined only if for each n > 1, the operator [(a + n)Id + A ] is invertible. The next 
step is then to search under what mathematical conditions the series defined by Eq. (A.5) converges. However, this 
step is quite complicated and therefore beyond the scope of this paper. 



A.2. Application 

In our case, we are only interested in obtaining the leading behaviour of our solutions near the surface. Therefore we 
will only solve the zeroth order equation (Eq. (A.4)) after having established the expressions for Yq and Ao. We start 
by defining x = 1 — C as the variable that will be used in the Frobenius series. The surface of the star then corresponds 
to x = and its interior to positive values of x. 

It is then necessary to choose a vector Y so that Eqs. (33)-(38) can be put in the form given by Eq. (A.l). This implies 
choosing the variables which are differentiated once with respect to x. Our choice is therefore Y — [II, vfi, Q = d x ^, \&]*. 
Having chosen the vector Y, it is then necessary to find the associated system of equations, by eliminating the variables 
(b, u e ,u^), and then to extract the zeroth order equation (see Eq. (A.4)). In fact, it is much simpler to do both steps 
simultaneously, given the complexity of Eqs. (33)- (38). 

Before giving the final result, it is important to point out that when N is not an integer, a mild singularity occurs 
on the surface of the star, due to the presence of fractional powers in the enthalpy, starting with x N+2 (Hunter, 2001). 
This in fact invalidates the use of Frobenius series in its present form from a strictly mathematical point of view, since 
these only use integer powers of x. This problem can be solved by including fractional powers in the solution, as is done 
in Christensen-Dalsgaard & Mullan (1994), the lowest one being x N+1 (this is not in contradiction with Christensen- 
Dalsgaard & Mullan (1994), for which the variable y 4 contains x N , because 7/4 includes d xx H in its expression whereas 
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our variables do not). As a result, the zeroth order equation remains unaffected and can therefore give the correct 
behaviour of the solution near the surface: 
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(A.6) 



This equation is based on the following development of the enthalpy near the surface: 
1 



H(x, 9) = H x {9)x + -H xx (9)x 2 + H N (0)x 



N+2 



0{X% 



(A.7) 



The characteristic equation is det(Ao — X.Id) = X 4 — NX 3 = 0. The eigenvalues are therefore a = —N and a = 0, 
the second value being triply degenerate. The first value is rejected because it leads to solutions that diverge on the 
surface of the star. The three remaining eigensolutions are bounded near the surface, which is in complete agreement 
with the results of Hurley et al. (1966), who applied the Frobenius method to the spherical case. By choosing a = 0, 
we also ensure that [(a + n)Id + Aq] is invcrtiblc for n > 1. These three bounded solutions and any of their linear 
combinations satisfy the following analytical constraint: 



An n = 



A(l -e)i?r ' 



(A.8) 



which is, in fact, equivalent to saying that Sp/po goes to zero on the outer boundary (where Sp represents the Lagrangian 
variation of the pressure). We can use the previous results to establish the behaviour of different quantities near the 
surface: 
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(A.9) 
(A.10) 
(A.11) 
(A.12) 
(A.13) 
(A.14) 



v is x 1 (this 



Eq. (A.8) shows that 8p/H N = o(l). The next relevant power in a power series expansion of Sp/H 
remains true even when N is not an integer since the first fractional power of Sp/H N is N + 1). By applying the 
equation Sp = c^Sp, it can also be shown that Sp/H 1 ^^ 1 = O(x). As a result we obtain the following behaviour for 
both Lagrangian perturbations: 



dp 
dp 



0(x N ) 



(A.15) 
(A.16) 



The results on p, p, 5p, 5p are interesting when we consider the equilibrium model. Since p <x H N and P cx H N+1 , 
we deduce that the leading behaviour of the equilibrium density and pressure are p Q = O (x N ) and P = O {x N+1 ), 
respectively. This implies that the ratio of the Eulerian density perturbation to the equilibrium density (p/p ) and the 
corresponding ratio for pressure both become unbounded as one approaches the surface of the star. This is problematic 
because the sum p a + A p cos(wi) (which corresponds to the total density) will periodically reach negative values close 
to the surface of the star for any non-zero amplitude A. However, the ratio of the Lagrangian density perturbation to 
the equilibrium density remains bounded as one approaches the surface, and the same applies to the pressure. This 
suggests that a Lagrangian description is physically more appropriate. 



Appendix B: Projection onto the spherical harmonic base 

B.l. Integral operators 

In order to project the fluid equations onto the harmonic basis, it is necessary to define a number of integral operators. 
The prototype to one of these operators is as follows: 



J%> (G) (0 = ff G(C, 6)doYF(0, <(>) {Y e m (9, dfi, 



(B.l) 
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where d£l = sm8d6d<fi, G is an arbitrary function, x* is the complex conjugate of x and J$, (.) is the operator. J$, (G) 
is a two-dimension array of indexes I and I' (the value of m is fixed) composed of functions depending on £ only. The 
remaining operators are given in the following table: 







{d e Y e m }* 


{D r™}* 


1 V 


Ju- (G) 
K%> (G) 


Jcf e , (G) 
L% (G) 
M%, (G) 


KcT t , (G) 
McT e , (G) 

m (G) 



If G is a real function than If v (G), J%, (G), Jc%, (G), Z$, (G) and N$, (G) are all real functions whereas K%, (G), 
(G), M^, (G), and McJJ, (G) are purely imaginary. There are symmetries between some of these operators: for 
example J%, (G*) = {Jc^ e (G)}*. The same applies for K%, (G) and Kcf t , (G), and for M£, (G) and Mc™ (G). 

In order to calculate these integrals, we use Gauss' quadrature. This gives accurate integrals when G is a "polyno- 
mial" of cos# (the coefficients of the polynomial can depend on (), for the operators 1$, (G), LjJ (G), M^, (G), 
Meg, (G) and N%, (G). For the operators (G), K$ (G), Jbg, (G) and Acg, (G), G needs to be of the form 
sin #P(cos 6>) where P is a polynomial. These integrals are calculated with i ros collocation points, where i ros is 
generally greater than L max , the harmonic resolution of the pulsations. 

Having defined the different integral operators, it is now possible to give explicitly the fluid equations projected 
onto the spherical harmonic basis. In what follows, we have used the following conventions: 



IZ,(G)ui = 1 £l%,(G)vf n , 



+NI 



\{G) = -I%, (G) + N%,(G). 



(B.2) 



It is also worth pointing out that in the following matrices, the summation on £' applies to an entire line of the matrix. 
For example, 



+7,7, (A) - J%, (B) ui 
-KrAG)+N%,(D) v e m 

is equivalent to: 

L L 

{ITAA)-JTAB)}u l m + {-K%(C) + N%,(D)}vi. 



(B.3) 
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B.2. Continuity equation 
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where we have made use of the following identities: 



(B.5) 



-l(£ + l)Y™ = d 2 99 Y™ + cot 6dgYl y 
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(B.6) 
(B.7) 



B.3. Adiabatic energy equation 
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B.4. Poisson's equation 



_i_ rm 



r 2 + rl 



[ -IS, (r 2 H N ^) 



CC m 
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(B.9) 



where we have made use of the Eq. (B.6). 



B.5. Euler's equations 
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Appendix C: The variational test 

The present formulation of the variational test is the same as that of Unno et al. (1989), apart from the following 
differences: we use the velocity rather than the displacement, hence the extra time derivatives; the star's volume is no 
longer spherical; the integral on the gravity wave energy is based on the effective gravity and uses the local vertical 
direction rather than e r ; the integral on the gravitational potential energy has been extended to infinite space. 

The different resultant integrals are given by the following explicit formulas and are calculated numerically using 
Gauss' quadrature in the angular direction and a spectral expansion in the radial direction (we use a radial resolution 
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of 101 points and an angular resolution of 200 points): 



/ Po\\v 

Jv 

I PoN 2 a \v 
Jv 

JV Po4 

[ Po Cl ■ (v* x v) dV 
Jv 

I 



■e g \ 2 dV 



dV 



I 

L 
I 



H 



N 



+ I«*l a £ + 2»{(«<)V} 



dV, 



NH 



N-l 



(N + 1)AH N ~ 1 



Ti 



J_ 

Tj r 4 r 2 

\u\ 2 dv, 



ku^dcH + u 6 d e H\ 2 dV 



2i n[ H N ( c _^l+ r ^l) 
Jv \ r C rr C J 



0_ + rgsm9 \ C 2 (u 6 r uj - ufufj + C 3 sin 9 (ujuj - u*uj) 



r 2 rQ 



|V*||W 



'V or V 2 



-f 

Jv 



r 2 +r 2 



7^H«c*l a + ilft*l a + — 



r 2 sin 9 2 



\d^\ 2 - ^n(d ( **d e t>)dv, 

r tq 



dV, 



(C.l) 
(C.2) 
(C.3) 

(C.4) 

(C.5) 



where dV = r 2 \r^\ s'm9d9d(d(j>, w£ — 3f v% = (vt ) etc. For the integral on the gravitational potential, it is 
useful to decompose infinite space into three domains: V U V2 U V3 = V^. V is the volume of the star, V2 is the 
volume comprised between the star and the sphere of radius 2, and V3 is the space outside the sphere of radius 2 (see 
Fig. 1). The integral on the first two domains is given by the expression above. For the third domain, it is based on the 
spherical harmonic decomposition of the gravitational potential. In empty space, a gravitational potential will take on 
the following form as it obeys the equation AW = and vanishes towards infinity: 



where the A 1 are constants. This form of W then leads to the following expression: 
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where r ext — 2 is the radius of the inner sphere of V3. This expression corresponds to the surface integral of Unno 
et al. (1989). 



